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ABSTRACT 

Computing  the  displacements  of  atoms  undergoing 
deformation  in  a  perfect  lattice  requires  the  use  of  the  so- 
called  Cauchy-Born  approximation.  Near  defects  such 
as  at  vacancies  or  dislocations,  this  approximation  does 
not  hold  and  a  computationally  costly  energy 
minimization  over  a  large  number  of  atoms  is 
unavoidable.  Presented  is  a  self-consistent  multiscale 
methodology  enabling  an  approximation  to  the 
displacements  near  defects  without  requiring  full  energy 
minimization  over  successive  load  increments.  It 
enables  description  of  the  fundamental  mechanical 
behavior  of  crystalline  materials  at  the  length  scale  of  a 
macroscopic  continuum  (i.e.,  millimeter  resolution) 
given  discrete  atomic  interactions  at  the  nanoscale  (i.e., 
angstrom  resolution).  The  basic  formulation  is  derived 
from  mathematical  homogenization  which  requires  the 
definition  of  a  representative  crystalline  volume  element 
containing  periodic  defects.  Numerical  simulations 
demonstrate  the  utility  of  the  framework  for  the 
particular  case  of  tungsten.  Elastic  stiffness  and 
energetic  properties  of  periodic  unit  cells  containing 
vacancies,  screw  dislocations,  and  low-angle  twist 
boundaries  are  computed  followed  by  demonstrative 
calculations  of  the  extension  to  atomistic  plastic  flow. 

1.  INTRODUCTION 

In  the  field  of  solid  mechanics,  one  of  the  key 
impetuses  for  multiscale  modeling  is  to  reduce  the 
empiricism  of  constitutive  models  to  an  irreducible  level. 
Typically,  this  is  the  level  of  atomistic  interactions.  The 
postulated  benefits  of  such  a  reduction  are  numerous. 
All  continuum  phenomenology  would  be  intrinsic  to  the 
atomistic  interaction  energy  either  in  terms  of  quasistatic 
atomic  forces  or  through  the  small  atomic  vibration 
periods  (femtoseconds).  Moreover,  no  assumption  of 
domain  simple  connectedness  would  be  required  so  that 
materials  may  develop  internal  interfaces  or  sliding 
boundaries  naturally. 

In  progressing  towards  these  benefits,  atomistic-to- 
continuum  scale  multiscale  techniques  continue  to  enjoy 
strong  consideration  for  modeling  nanoscopic  structures 
that  are  difficult  to  treat  with  atomistics  alone  (Tadmor  et 
al.,  1996;  Rudd  &  Broughton,  2001;  Ghoniem  et  al., 
2003).  With  the  ever  increasing  power  of 
microprocessors,  one  might  eschew  multiscale  methods 
and  instead  pursue  brute  force  large-scale  atomistic 
calculations  (cf.  Seppala  et  al.,  2004)  to  evoke  continuum 


behavior.  However,  at  the  current  rate  of  increase  in 
performance,  many  years  still  remain  before  computers 
will  be  capable  of  simulating  the  response  of  enough 
atoms,  over  a  time  scale  of  sufficient  duration  to 
realistically  address  design  of  “engineer- able”  materials. 

Despite  the  successes  of  atomistic-to-continuum 
modeling  methods,  significant  limitations  still  remain. 
Most  notable  among  these  is  the  requirement  of 
kinematic  approximations  on  deforming  crystals  that 
preclude  multiscale  methods  from  handling  important 
solid  mechanics  problems  in  general  electromagnetic 
materials  (Cole,  et  al.,  2003).  The  kinematic 
approximation,  also  called  the  Cauchy-Born  rule  (cf. 
Ericksen,  1984),  reduces  the  number  of  computational 
degrees  of  freedom  by  approximating  the  displacements 
of  lattice  points  through  affine  maps  using  the 
deformation  gradient.  As  a  mathematical  tool  for 
multiscale  modeling,  the  Cauchy-Born  rule  is 
indispensable  as  it  unifies  the  behavior  of  all  atoms  in  a 
perfect  crystal,  effectively  enabling  the  replacement  of 
many  atom  degrees  of  freedom  with  a  single  second  rank 
tensor.  Crystals  that  are  amenable  to  the  Cauchy-Born 
approximation  are  typically  monocrystalline,  covalent,  or 
centro  symmetric . 

For  imperfect  crystals,  or  more  specifically  crystals 
that  lack  an  axis  of  inversion  symmetry  (Ericksen,  1984), 
no  similar  approximation  currently  exists.  These  types 
of  crystals  usually  are  complex,  ionic,  or  contain  lattice 
defects.  Most  engineered  materials,  particularly  of 
interest  to  the  Army  and  long-term  commercial 
application,  contain  critical  components  whose  crystals 
are  of  this  type. 

In  this  paper,  we  develop  a  numerical  method  to 
enable  multiscale  modeling  of  imperfect  crystals.  The 
method  is  based  on  but  computationally  less  expensive 
than  an  asymptotic  homogenization  technique  developed 
earlier  by  the  authors  for  hyperelastic  semiconductor 
crystals  (Chung  &  Namburu,  2003;  Chung,  2004).  We 
demonstrate  its  use  for  zero  temperature  mechanics  of  a 
crystalline  transition  metal  solid  with  comparisons  to 
direct  molecular  static  results.  Under  strain,  we  examine 
point,  line  and  plane  defects.  Additionally,  we  present 
initial  results  for  the  extension  of  the  method  to  handle 
bulk  plasticity  using  only  the  atomistic  information  for 
the  constitutive  properties  through  the  introduction  of  an 
evolving  reference  configuration. 

The  solid  mechanical  properties  of  interest  are  the 
effective  elastic  stiffness  and  strain  energy  density. 
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Section  2  overviews  the  multiscale  homogenization 
equations  applicable  to  nonlinear  hyperelastic  problems. 
Section  3  introduces  the  notion  of  discrete  atoms  in  the 
context  of  homogenization.  This  approach  is  extended  to 
numerous  periodic  defect  examples  in  Section  4, 
followed  by  discussions  and  conclusions  in  Section  5. 


2.  MATHEMATICAL  DEVELOPMENTS 


2.1  Governing  equations 

Two-scale  homogenization  involves  the  self- 
consistent  solution  of  two  coupled  equations  derived 
from  the  statement  of  conservation  of  linear  momentum. 
In  brief,  for  a  given  energy  density  F  in  reference 

configuration  B 0 ,  current  configuration  B  at  time  t , 
with  reference  macroscopic  coordinates  X  and  current 
coordinates  xa  =  xa  (XA ,  t ) ,  the  respective  microscopic 

coordinates  Y  and  ya  =  ya  (YA ,  t) ,  coarse-  and  fine-scale 
displacements  u  and  v,  the  corresponding  macroscopic 
and  microscopic  deformation  gradients 

F]  -du  /  dXA  +  5\  and  f]  -  dya  /  dYA  +  Sa. ,  and  the 
assumption  of  additive  decomposition  of  displacements 
u“  =ua+u  =  ua  +  sva  (Chung,  2004),  we  have 


1  <•  j -d'F  d(Su °) 
Y  1 1  dF"B  dXB 


dVdY  =  |  T'gJm'dA 

dV 


+  \B°gab5ubdV,  (vSua) 


(1) 


1  j-  j -d'F  d(8v) 
Y  1 1  dF“  8Yb 


dVdY  =  0. 


(v£v“ ) 


(2) 


du a  dv° 

F  = - + - +  8 

dXA  dYA 


fVu  A  ' 

— st  +  s: 

KdXA  a  ‘y 


Fa  =  FlFa  , 


(4) 


where  Ft.  is  the  deformation  gradient  under 
microscopically  homogeneous  conditions  and 
F A  =  F~laa FaA  depends  upon  the  gradient  of  v  and  thus 
accounts  for  heterogeneity  due  to  defect  fields.  Barred 
Greek  indices  denote  a  third  configuration  B  for  the 
material  associated  with  the  multiplicative 
decomposition  in  the  last  of  (4). 


2.2  Homogenization  with  an  evolving  intermediate 
frame 

Consider  now  a  continuum  elastic-plastic  material 
where  the  deformation  gradient  F  is  decomposed  into  a 

lattice  part,  FL ,  and  a  plastic  part,  FP  (Bilby  et  al., 
1957;  Teodosiu,  1969): 

F  =  FlFp  ,  F  =  FLaFPaA .  (5) 

This  presumes  the  existence  of  an  intermediate 
configuration  B  ,  with  FP  the  tangent  mapping  between 
Bq  and  B  ,  and  FL  the  tangent  mapping  from  B  to  B . 

The  plastic  map  FP  reflects  contributions  to  the  total 
deformation  that  leave  the  lattice  unperturbed,  such  as 

the  flux  of  mobile  dislocations.  The  lattice  map  FL 
includes  all  other  kinematic  contributions,  including  rigid 
body  motions,  elastic  fields,  and  residual  elastic  fields 
attributed  to  the  lattice  defects. 


2.3  Plasticity 


This  set  of  equations  defines  conventional  quasistatic 
asymptotic  homogenization.  The  coarse-scale  Eq.  (1) 
and  fine-scale  Eq.  (2)  are  coupled  through  the 

constitutive  dependency  F  =  *F  (F  (du  /  dX,  5v  /  9Y))  . 


In  an  elastic-plastic  body,  it  is  typically  assumed  that 
the  free  energy  depends  only  upon  that  part  of  the 
deformation  that  affects  the  lattice,  in  other  words, 

F  =  JP[F(  F,  FP  X0)  =  ^(FL  X0)9  (6) 


In  equations  (1)  and  (2),  ua  represents  the 
displacement  that  would  exist  in  a  microscopically 

homogeneous  medium  and  ua  accounts  for  fine-scale 
heterogeneity,  with  corresponding  fine-scale 

representation  va .  The  corresponding  microscopic 
decomposition  is 

a  -1  a  — a  .  ~a  (  -t-ia  c*a  \  t tA  ,  ~a 

V  =£  U  =v  +  v  =[Fa-Sa)Y  +  v  ,  (3) 

with  v a  the  microscopic  displacement  arising  from  the 
projection  to  the  fine-scale  of  the  macroscopic 

deformation  gradient  FA  .  From  here  we  may  write  the 
decomposition 


where  F  is  measured  per  unit  intermediate  volume  on 
B ,  and  where  we  have  also  included  the  dependence  of 
the  free  energy  on  absolute  temperature  6 . 

Additionally,  %  is  a  vector  of  internal  state  variables 
accounting  for  deviations  in  stored  energy  from  that  of  a 
perfect  lattice.  Notice  that  generally,  FL  and  FP  of  Eq. 

(5)  and  F  and  F  of  Eq.  (4)  are  four  distinct  deformation 
mappings,  with  corresponding  configurations  of  the  body 
depicted  in  Fig.  1(a).  The  former  two  denote  lattice  and 
plastic  deformations;  the  latter  two  represent  micro- 
homogeneous  and  micro-heterogeneous  deformations, 
respectively. 

3.  ATOMIC  HOMOGENIZATION 
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3.1  Basic  Considerations 


The  equations  thus  far  have  referred  to  continuum 
conservation.  To  introduce  the  notion  of  discrete  atoms, 

define  spatial  positions  of  atoms  in  configuration  B 

in  Cartesian  coordinates: 


z{j) 


r%a  r-w  a  a 

8 aZ(j)  +<!(])■■ 


(7) 


with  a  displacement  vector  between  intermediate 
and  current  states.  Let  and  denote  interatom 

vectors  in  respective  configurations  B  and  B  ,  i.e., 

*C*>  =  %> _  h)  ’  V>  =  z« _  zo>  • 

Making  contact  with  the  homogenization  theory  of 
Section  2,  we  define 


a  _  t-t  L  a  ry  cl  ~a 

Z(j)  0k)a  ( k >  V(j)  ' 


(9) 


The  first  term  in  Eq.  (9),  F^aZa^ ,  accounts  for  the 

uniform  projection  over  each  periodic  cell  of  the 
macroscopic  lattice  deformation  field  to  the  fine  scale 

(i.e.,  the  Cauchy-Born  approximation),  and  is  the 

discrete  atomistic  analog  of  the  perturbation  in 
displacement  due  to  microscopic  heterogeneity. 


I  B  i 


Figure  1.  Deformation  mappings  and  configurations  (a), 
and  purely  plastic  deformation  (b). 

After  considerable  equation  derivation,  we  arrive  at  the 
analogous  atomistic  set  of  homogenization  equations 
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\rgabSubdA  +  \BagabSuhdV- 
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The  mapping  FP  represents  the  cumulative 
deformation  of  the  material,  but  not  the  lattice,  due  to 
dislocation  glide  through  the  unit  cell.  For  a  fixed 
control  volume,  once  atoms  have  convected  through  the 

lattice  due  to  FP  ,  BQ  and  B  appear  identical  in  the  fine 

scale  domain  if  no  dislocations  are  created  or  destroyed 
within  the  volume  element,  yet  net  deformation  of  the 
material  will  have  taken  place  at  the  coarse  scale  (Asaro, 

1983).  For  a  purely  plastic  process  with  FL  =  1,  new 
atoms  would  enter  the  control  volume  (i.e.,  the  atomic 
unit  cell)  in  the  identical  locations  as  the  old,  to  replace 
those  that  exited  due  to  plastic  flow.  The  concept  is 
illustrated  in  Fig.  1(b). 

3.2  Elastoplastic  Atomistic  Homogenization 

We  assume  a  free  energy  potential  depending  only 
upon  the  relative  positions  of  all  atoms  (i.e.,  lattice 

statics)  Y  =  W  (q0> ,  Z(y> )  =  W  (r(lu) ,  r<IU> ,  r<M> ,  ...r(N_nN> )  . 

Under  isothermal  conditions,  we  may  rewrite  this  as 

=  (10) 

where  %  accounts  for  effects  on  stored  energy  due  to 
deviations  from  homogeneous  elasticity  at  the  fine  scale. 


1 

Y 


I  JA>vo> 


d(Sua ) 
dXB 


dVdY. 


where  Ku  =  -tfr/dq^dF^  ,  =  d2F/dq{k)dqb{l) , 

c:=J~F»'c:f™  and  =Dlj)ba=r'F-pi^.  In 

subsequent  calculations,  (11)  is  solved  for  the  (atomic) 
inner  displacements  and  (12)  can  be  solved  using  a 

standard  numerical  procedure  such  as  a  finite  element 
method. 


3.3  Inelastic  fields  and  defect  kinetics 

For  simplicity,  temperatures  are  presently  omitted. 
Therefore,  additional  kinetic  equations  are  required  to 

advance  the  plastic  deformation  gradient  FP  and  internal 
variables  In  the  continuum  description,  these  are 
written  as 

FP  =  FP  (f,Fp,|,$),  (13) 

I  =  |(f,Fp,|,6»)  .  (14) 
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One  can  formulate  an  approximate  expression  for  the 

rate  of  the  mapping  FP  in  terms  of  the  dislocation 
velocity,  line  length,  and  orientation  (Teodosiu,  1969). 


A  similar  approach  has  been  used  by  Zbib  and  co¬ 
workers  in  the  context  of  discrete  dislocation  plasticity 

(Zbib  &  de  la  Rubia,  2002).  The  determination  of  FP 
is  nontrivial.  In  the  absence  of  temperatures,  we  must 
assume  specific  forms  in  section  4.3. 

3.4.  Numerical  Implementation 

In  cases  of  monotonically  loaded  crystals,  we 
employ  the  improved  computer  procedure  of  (Clayton  & 
Chung,  2006)  wherein  a  full  solution  to  Eq.  (11)  is 
evaluated  only  for  the  first  load  step  and  subsequent 
steps  are  obtained  via  the  minimization  of  a  scalar 

numerical  parameter  r/1 ,  through  a  line  search  algorithm. 
This  is  achieved  by  a  modification  to  (9)  such  that  it  now 
reads 


Z(J)  -  F(jk)ct  Z{k)  +  V(j) 


(15) 


where  the  index  i  is  the  load  increment  and  v^'  is 
defined  by  v^1  =  rj 'v^  . 


4.  APPLICATION:  DEFECTS  IN  TUNGSTEN 

The  present  formulation  is  applied  to  study  the 
mechanical  behavior  of  tungsten  (W),  a  BCC  transition 
metal  of  relatively  high  mass  density.  Its  combination  of 
high  density  and  strength  make  it  an  attractive  material 
for  use  in  defense  applications  such  as  ordnance  (Zhou  & 
Clifton,  1997;  Clayton,  2005).  Specifically  computed 
here  are  the  nonlinear  elastic  responses  of  deforming  W 
crystals  containing  (i)  vacancies,  (ii)  screw  dislocations 
of  like  sign  and  screw  dislocation  dipoles,  or  (iii)  low- 
angle  twist  boundaries.  In  section  4.2,  unit  cells  are 
deformed  in  uniaxial  stretch  to  2.5%  elongation  or  larger. 
The  computational  method  is  demonstrated  by  direct 
comparison  with  conjugate  gradient-based  lattice  statics. 
This  is  followed  by  results  in  section  4.3  of  larger- 
deformation  shear  simulations  (to  10%  applied  shear 
strain)  reported  for  the  case  of  screw  dislocation  glide  in 
BCC  W,  in  which  the  evolving  intermediate 
configuration  is  updated  assuming  monotonic  single  slip 
occurs  as  resolved  shear  stresses  exceed  the  Peierls 
threshold  (cf.  Hirth  &  Lothe,  1982). 

4.1  Atomistic  potential:  tungsten 

Without  loss  of  generality  in  the  underlying  method,  we 
employ  an  empirical  N  -body  potential  specifically 
developed  for  transition  metals  (Finnis  &  Sinclair,  1984) 

in  order  to  estimate  the  free  energy  potential  W  of  (10). 
Duesberry  &  Vitek  (1998)  exercised  the  Finnis- Sinclair 
potential  to  model  screw  dislocation  core  structures, 
construct  generalized  stacking  fault  energy  surfaces,  and 
predict  plastic  slip  anisotropy  (i.e.  tension-compression 


asymmetry  in  the  yield  function)  in  multiple  BCC 
transition  metals. 

4.2  Static  defects  in  tungsten 

The  calculations  are  performed  to  determine  the 
effective  tangent  stiffness  of  a  microscopic  (i.e., 
atomistic)  unit  cell.  A  sample  atomistic  unit  cell 
representative  of  all  calculations  described  in  the  present 
work  is  shown  in  Fig.  2.  The  unit  cell  dimensions 

areL  xL  xL  =  a'\J3N  xa^j6N  xa\2N  ,  where  N., 

N2,  and  N3  are  respectively  the  number  of  repeating 

planes  stacked  in  the  [111]-,  [112]-,  and  [In¬ 
directions.  Periodic  boundaries  in  the  form  of  wrap¬ 
around  conditions  are  applied  along  all  faces  of  the  unit 
cell  in  the  usual  manner  (cf.  Vitek,  1976). 


Figure  2.  Atomistic  scale  unit  cell  for  perfect  BCC 
lattice. 

First,  we  examine  the  response  of  W  containing  unit 
cells  with  immobile  but  deformable  lattice  defects,  i.e. 

the  local  deformation  gradient  F  =  FL  =  FF  in  the 
context  of  Fig.  1.  Fundamental  material  effects  in  the 
context  of  lattice  statics  calculations  have  been  examined 
by  (Vitek,  1976;  Duesberry  &  Vitek,  1998). 

The  initial  atomic  coordinates  are  found  using  a  two- 
step  procedure:  first  the  linear-elastic  solution  for 
displacement  field  of  the  defect  is  applied  to  the  atoms, 
then  a  conjugate  gradient  algorithm  (Plimpton  & 
Hendrickson,  1993)  is  used  to  find  local  minimum 
energy  atomic  positions.  The  homogenization  method 
then  proceeds  to  find  the  changes  in  local  states  with 
applied  deformation.  The  applied  (i.e.,  coarse-scale) 
deformation  gradient  field  (in  conjunction  with  fine-scale 
periodicity)  is  uniaxial  stretching  in  the  [1 1 1] -direction 

over  a  range  of  1.000  <  Fn  <  1.025 .  We  verify  the 

computed  configurations  by  comparing  to  incremental 
conjugate  gradient  (CG)  energy  minimization  used  in 
conventional  lattice  statics. 

4.2.1  Point  Vacancy 

The  initial  configuration  for  the  point  vacancy  is 
constructed  by  removing  one  atom  closest  to  the  centroid 
of  the  unit  cell.  The  defect  density  in  this  case  is  defined 

as  the  volume  fraction  of  missing  atoms,  i.e.,  pd  =  1  /  N  , 
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0.0020 


(d) 

Figure  3.  Vacancy  configuration:  (a)  strain  energy 
comparison  between  computational  methods,  (b)  defect 
energy,  (c)  elastic  stiffness,  and  (d)  Zener  anisotropy. 


where  N  is  the  total  number  of  atoms  prior  to  vacancy 
creation  ranging  from  2016  to  32256.  Fig.  3(a)  shows 
the  defect  energy  comparison  with  CG.  The  agreement 
is  excellent.  In  Fig.  3(b-d)  are,  respectively,  the  defect 
energy  under  variations  of  the  defect  density,  the  elastic 


stiffness  C1111,  and  the  Zener  anisotropy  factor  A,  all 
under  parametric  variations  of  the  defect  density.  The 
Zener  anisotropy  factor  A  is  defined  by  (cf.  Hirth  & 
Lothe,  1982) 


Note  that  defect- free  W  is  nominally  isotropic  at  null 
deformation,  i.e.,  ^4  =  1.00,  and  that  the  anisotropy 
increases  dramatically  as  the  lattice  is  stretched.  First- 
principles  (i.e.,  quantum-mechanical)  calculations  have 
indicated  a  departure  from  isotropy  in  W  at  large 
pressures  (Ruoff  et  al.,  1998). 


4.2.2  Screw  Dislocation 

The  dislocation  tangent  line  and  burgers  vector  b  for 
the  screw  dislocation  example  are  oriented  along  the 
[1  Indirection  and  pass  through  the  centroid  of  the  unit 

cell,  the  latter  having  a  magnitude  of  b  =  |b|  =  yjla/ 2  . 
The  scalar  dislocation  density  is  defined  as  the  defect 


line  length  per  unit  reference  volume,  i.e., 
pd  =  1  /  (Z0Z3 )  •  The  defect  energy  per  atom,  as  shown 

in  Fig.  4(a),  is  computed  accurately  by  our 
homogenization  scheme  as  verified  by  the  incremental 
CG  solutions.  We  see  a  linear  increase  in  stored  defect 

energy  with  Fn ,  and  a  roughly  linear  increase  with 

increasing  pd .  Gibeling  &  Nix  (1980)  discussed,  from 
the  standpoint  of  discrete  dislocation  modeling,  how  the 
strain  energy  supported  by  dislocations  may  be  amplified 
by  applied  external  deformations,  and  Clayton  (2005) 
assumed  a  linear  dependence  of  stored  elastic  energy  on 
dislocation  density  in  a  continuum  crystal  plasticity 
model  of  single  crystalline  W.  From  Fig.  4(c)  and  Fig. 
4(d),  the  dislocation  density  tends  to  accelerate  the 

decrement  in  C1111  with  increasing  applied  stretch, 
whereas  A  tends  to  be  suppressed  by  increasing 
dislocation  content  as  the  stretch  increases.  Although  the 
trend  of  decreasing  elastic  modulus  with  dislocation 
content  has  been  reported  from  physical  experiments 
(Smith,  1953)  and  analytical  continuum  modeling 
(Lebedev,  1996),  it  has  not  been  emphasized  previously 
in  continuum  plasticity  models. 


4.2.3  Screw  Dislocation  Dipole 


The  dislocation  tangent  lines  for  the  screw 
dislocation  dipole  are  oriented  along  the  [1 1  Indirection. 
The  pair  with  opposite  burgers  vector  signs  are  located  at 

(/L2,/4L3)  and  (%L2,%L3) ,  maintaining  a  minimum 

separation  distance  of  (X^vXA)  between  all 
dislocations  whether  within  the  cell  or  in  images. 


The  dislocation  density  pd  =  2  /  (C2L3 ) ,  and  the  net 

Nye  (1953)  tensor  vanishes,  as  the  burgers  vectors  cancel 
out,  meaning  the  dislocation  density  is  ‘statistically 
stored’  in  the  sense  of  Ashby  (1970).  The  defect  energy 
computed  with  our  homogenization  scheme  is  plotted  in 
Fig.  5(a);  again,  in  good  agreement  with  CG  simulations. 

From  Fig.  5(b),  we  see  that  Ed  increases  linearly  with 
pd ,  and  it  increases  linearly  with  Fu  .  As  shown  in  Fig. 

5(c),  the  stiffness  coefficient  C1111  decreases  with 
increasing  defect  density  and  increasing  stretch. 
Anisotropy  A  (Fig.  5(d))  increases  with  stretch  and 

increases  with  dislocation  content  at  low  values  of  pd 

(<  0.013  /nm2 ),  but  decreases  as  pd  is  increased 
further. 


4.2.4  Low  Angle  Twist  Grain  Boundar 


The  final  defect  configuration  examined  is  a  low 
angle  twist  grain  boundary.  We  describe  this  boundary 
type  using  the  disclination  concept  (Li,  1972;  Clayton  et 
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— • -  Homogenization,  p  =  2.6(1 0)'2/nm2 

O  Conjugate  gradient,  p  =  2.6(10)'2/nm2 


(a) 


Figure  4.  Screw  dislocation:  (a)  strain  energy 
comparison  between  computational  methods,  (b)  defect 
energy,  (c)  elastic  stiffness,  and  (d)  Zener  anisotropy  (e). 


defect  density  pd  is  increased  from  0  to  0.026/nm2.  At 
an  applied  stretch  of  Fn  =  1.025 ,  C1111  increases  over  the 

same  range  of  pd  from  414  GPa  to  419  GPa. 
Anisotropy,  shown  in  Fig.  6(d),  is  suppressed  slightly 
with  increasing  disclination  content. 

4.3  Multiscale  elastoplasticity 

Since  the  flux  of  dislocations  leaves  the  crystal 
lattice  relatively  unperturbed  upon  complete  glide  across 
the  periodic  unit  cell,  we  appeal  to  the  notion  of  an 
evolving  relaxed  intermediate  or  natural  configuration 
(Bilby  et  al.,  1957;  Teodosiu,  1969;  Asaro,  1983)  and 
introduce  it  explicitly  into  the  formulation  in  this  paper 
for  plasticity.  In  the  approach,  finite  inelastic 
deformation  is  accounted  for  by  the  evolution  of  this 
relaxed  configuration,  as  dictated  by  the  dislocation  flux 
through  the  volume  element.  The  instantaneous 
thermoelastic  response  of  the  element  is  measured  with 
respect  to  the  natural  configuration,  and  driving  forces 
for  defect  propagation  are  calculated  at  the  current  state 
or  time  instant,  in  order  to  project  the  intermediate 
configuration  forward  temporally. 

The  full  decomposition  F  =  FLFP  =  FF  of  Fig.  2 
applies.  Initial  conditions  for  the  unit  cell  are  identical  to 
those  of  Fig.  4(a).  We  apply  pure  shear  deformation  in 
the  1-3  plane,  over  the  range  0.00  <  Fu  <0.10.  The 

multiplicative  decomposition  F  =  FLFP  yields 


al.,  2006).  The  defect  misorientation  across  the  (1 10)- 
plane  is  co-  0.2  radians ,  a  value  that  ensures  a  stable 
local  configuration  for  the  initial  atomic  arrangement. 
The  Frank  vector  is  oriented  parallel  to  [1  10] ,  the 
disclination  line  i  is  oriented  along  [111],  and  the  axis  of 
twist  passes  through  the  cell  centroid.  The  disclination 

line  density  is  pd  =  1  /  (L2L3 ) .  The  defect  energy  results 

under  deformation  are  again  validated  via  comparison 
with  an  incremental  conjugate  gradient  solution  in  Fig. 
6(a).  The  agreement  is  once  again  very  good.  From  Fig. 

6(b),  Ed  increases  nonlinearly  with  disclination  line 
density;  that  is,  it  increases  nonlinearly  with  grain 
boundary  area  per  unit  volume. 

Note  that  we  are  essentially  simulating  a  periodic 
array  of  such  low-angle  grain  boundaries,  as  have  been 
known  to  appear  in  pure  W  subjected  to  severe  plastic 
deformation  (Valiev  et  al.,  2002;  Wei  et  al.,  2005).  In 
contrast  to  the  other  classes  of  defects  examined  in  the 
present  work,  the  defect  energy  decreases  with  applied 

stretch,  and  the  stiffness  component  Cim  increases  with 
increasing  defect  density.  From  Fig.  6(c),  at  null  applied 

strain,  C1111  increases  from  514  GPa  to  520  GPa  as  the 


Fi3=r  =  rL+rP,  0?) 

where  yL  is  the  lattice  shear  associated  with  the  external 

stress  and  yP  is  the  cumulative  plastic  slip,  assumed  in 

our  idealized  problem  here  to  occur  as  a  result  of 
dislocation  glide  on  (110)-  planes.  Plastic  slip  is 
assumed  to  follow  the  simple  kinetic  relation 


,  ,  y[ 

(2  7tlab)x  j 

r  =r~r„  - — ii  +  sm 

2  l 

V 

_  (r-rl-ab/4)  JJ 

for  (\/y  >  y*  j,  and  yP  =  0  otherwise,  where  K  is 
the  initial  yield  strain,  y[  is  an  additional  lattice  strain 

required  to  overcome  the  Peierls  barrier,  and  a  -\j L3 

scales  the  amount  of  shear  strain  accumulated  in  the  unit 
cell  due  to  the  slip  of  a  dislocation  over  the  distance  of 
one  burgers  vector.  The  sinusoidal  character  of  (20)  is  in 
accordance  with  the  Peierls  model  of  periodic  lattice 
resistance  stress  discussed  by  Hirth  &  Lothe  (1982).  The 
simulations  discussed  here  assume  that  the  primary 
dislocation  remains  fixed  in  the  unit  cell  (in  the  atomistic 
domain),  and  that  no  defect  generation  or  annihilation 
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occurs,  as  other  screw  dislocations  convect  through  the 
lattice  and  cause  an  increase  in  plastic  shear  yP . 


Figure  5.  Screw  dislocation  dipole:  (a)  strain 
energy  comparison  between  computational  methods,  (b) 
defect  energy,  (c)  elastic  stiffness,  and  (d)  Zener 
anisotropy. 


Following  the  molecular  statics  calculations  of  Vitek 
(1976),  we  choose  yy  =  0.020  and  y[  =0.009,  the 

former  a  shear  strain  required  to  cause  [111](1  10)- 
screw  dislocation  core  shuffling  in  BCC  W,  the  latter 
following  his  result  that  a  total  shear  strain  of  0.029  is 
required  to  sustain  steady  dislocation  glide. 


Upon  initial  yield,  the  stress-strain  curves  follow  Peierls- 
type  oscillations  (a  consequence  of  using  (18)),  with  each 
period  corresponding  to  the  passage  of  a  single  [111]- 
screw  dislocation  across  the  unit  cell,  gliding  on  a  (110)- 

plane.  As  the  dislocation  density  pd  is  inversely  related 
to  the  dimension  L3  of  the  unit  cell,  the  larger  the 
dislocation  density,  the  smaller  the  unit  cell  size,  the 
larger  increment  in  shear  strain  yP  associated  with  the 
passage  of  each  dislocation,  and  the  longer  the  period  of 
shear  stress  oscillations  in  Fig.  7(a).  In  the  absence  of 
such  oscillations,  the  material  behavior  would  be  elastic- 
perfectly-plastic.  The  total  energy  per  atom, 

Ef  =  E  /  N  -Uc ,  is  shown  in  Fig.  7(c).  At  zero  applied 
shear,  differences  among  energy  curves  are  due  purely  to 
differences  in  defect  concentration,  while  the  energy 

oscillations  at  larger  shear  deformations,  Fu  >  0.020 ,  are 

associated  with  the  lattice  deformation  fields  (FL) 
appearing  in  conjunction  with  the  applied  stress. 

5.  DISCUSSION  &  CONCLUSIONS 

A  multiscale  method  predicated  on  asymptotic 
homogenization  has  been  developed  and  implemented. 
The  technique  is  based  on  conservation  principles 
formulated  using  deformation  gradients  at  both  micro 
and  macro  scales.  Aside  from  self-periodicity  near  the 
boundaries,  no  restriction  is  placed  on  the  atomic 
structure  or  for  continuity  with  the  supporting 
continuum.  Thus  the  formulation  directly  accounts  for 
the  effects  of  lattice  defects  on  the  homogenized 
mechanical  properties,  as  well  as  the  lattice-preserving 
kinematics  of  finite  plastic  deformation  associated  with 
dislocation  flux,  all  as  functions  of  the  state  of 
deformation.  From  a  computational  standpoint,  the 
homogenization  method  has  been  demonstrated  to 
accurately  reproduce  lattice  statics  in  terms  of  prediction 
of  minimum  energy  in  statically-deforming  lattices 
containing  defects.  It  enables  the  consideration  of  very 
large  continuum  scales  with  finite  nontrivial  atomic 
structures  embedded  within. 


The  1-3  shear  component  of  the  first  Piola-Kirchhoff 

stress,  Pu ,  is  shown  in  Fig.  7(a),  where  components  of 
the  stress  tensor  P  are  computed  from 

dE 

P  A  =JP  l - FP~lA  .  (19) 

dFLa ' 

a 

Prior  to  initial  yield,  slight  differences  in  slope  of  the 
stress-strain  response  are  evident  among  the  three  curves, 
each  curve  corresponding  to  a  different  dislocation 

density  pd .  Such  differences  are  due  in  part  to  the 
decrease  in  effective  elastic  shear  modulus,  C1313  =  C"  , 
with  increasing  defect  density,  as  is  seen  in  Fig.  7(b). 


The  homogenization  method  was  implemented 
numerically  to  study  the  nonlinear  elastic  response  of 
BCC  tungsten  lattices  containing  periodically-distributed 
vacancies,  screw  dislocations  and  dipoles,  and  low-angle 
twist  boundaries.  The  material  results  showed  that 
general  trends  of  defect  structures  subject  to  strain  are  in 
agreement  with  past  observations  and  with  strong 
quantitative  agreement  with  computed  fundamental 
material  properties.  Both  static  and  initial  dynamics 
results  show  conclusive  agreement  with  alternative,  more 
costly  approaches. 
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Figure  6.  Twist  disclination:  (a)  strain  energy 
comparison  between  computational  methods,  (b)  defect 
energy,  (c)  elastic  stiffness,  and  (d)  Zener  anisotropy. 
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Figure  7.  Elastic-plastic 
homogenization:  shear  stress 

(a) ,  elastic  shear  stiffness 

(b) ,  and  total  relative  elastic 
energy  (c). 
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A  Multiscale  Modeling  Method  for  Deformations  on 
Atomic  Lattice  Defects  and  Application  to  Plasticity 


Overview 


1.  Background:  crystal  kinematics 

2.  Homogenization 

Length  scales  and  coordinate  systems 

Elastic  continua 

Elastic  -  plastic  continua 

Atoms  -  discreteness  at  the  fine  scale 

3.  Defects  in  Body-Centered-Cubic  tungsten 

Atomistic  potential 

Numerical  results:  vacancy,  dislocation,  twist  boundary 
Numerical  results:  screw  dislocation  with  continuum  plasticity 


Deformation  kinematics  of  crystalline  solids 


Perfect  lattice  [Cauchy,  1828;  Born,  1914;  Ericksen,  1984;  Zanzotto,  1992] 


Continuum  mechanics:  dx=FdX 
Born-rule  mechanics:  da=FdA 


Imperfect  lattice:  Born  rule  does  not  apply  ,  _  ^  ,  ,  ,  _  fm 

K  y  c/x=Fc/X  but  da  t  FdA 


oooooo 

OOOQOO 

oooooo 

oooooo 

oooooo 

OOOOO^J^ 


dX  "  dx 


Multiscale  scheme: 
coordinates  and  configurations 


Externally  unstrained  configuration 


Deformed  configuration 


[Bensoussan  et  al.,  1978;  Sanchez-Palencia,  1980;  Bakhvalov  &  Panasenko,  1989;  Chung  &  Namburu,  2003;  Chung,  2004] 


Asymptotic  homogenization 
with  geometrically-nonlinear  elasticity 


Kinematics  (displacements) 


Hyperelasticity 


*F  =  !F(F) 


Substituting,  integrating  over  Y,  and  lim  s  —*■  0  gives 


dV  8(Sv“) 


YJyJrdF"B  dY‘ 


dVdY  =  0  (v<Jv") 


coupled  micro-  and 
macro-  equations; 
solve  concurrently 
for  va  and  ua 


1 


Incorporation  of  finite  plasticity  at  coarse  scale  - 


Kinematics  (displacements) 


Deformation  gradient 


Fine  and  coarse  scale  equilibrium  from  virtual  work  principle 


lim^  ->  0  I 
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Incorporation  of  atomic  kinematics  at  fine  scale  - 


Discrete  kinematics  (Cartesian  space) 


Fine  scale  equilibrium  (after  some  algebra) 


C 'a  rya  . 

:  oZ,  v  +q, 

a  (m)  J-h 


|R/  w\  —  Z/,\  —  Z/  v 

w  (a) 


ir(a\b)~Z(b)  Z(a) 


atom 

displacements 


relative  atom 
positions 


\z:  N  =F,  \  Z, 

(m)  {mnja  {n 


Coarse  scale  equilibrium 


lattice  mapping 


I V)  =  ¥L^(a\b)  +  V{6)  -  V(a)  =  FLR  +  f 


New  computational  method 
for  global  energy  minimization 


•  Incremental,  iterative  quasi-static  formulation  (longer  time  scales) 


•  Reduce  full  relaxation  problem  of  atom  DOFs  to  single  parameter  line  search 
with  bisection  (predict  structure  of  defect  cores  without  molecular  dynamics) 


•  Comparable  overall  cost  (DOF,  solution  time)  to  lattice  statics 


Numerical  examples 


Tungsten  with  Embedded  Atom  Method:  unit  cell  calculations  for  effective  properties 

-  Point  vacancy 

-  Screw  dislocation 

-  Low  angle  grain  boundary  (twist  disclination) 

-  Screw  dislocation  with  plastic  shear  kinematics 


Multiscale  Concept 


BCC  Unit  Cell 


Implementation  in  homogenization  scheme  (analytical  derivatives  and  chain  rule) 


Analytical  derivatives  ■  TTb 

substituted  into  equilibrium  M^(%b  ~ 1 =  ^(jk)ab^(k 

equation: 


Potential  used  by  Vitek  [1998]  and  others  for  dislocation  core  structures  in  W, 
but  has  drawbacks  for  atoms  too  tightly  packed  [Ackland  &  Thetford,  Phil.  Mag.  A,  1987]. 


Vacancy  -  [100]  tension 
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Numerical  Validation 


Defect  Configuration 


X  (110) 


Defect  Energy 


0.0020 


0.0015 


E 

o 

ro 

I 


0.0010  - 


0.0005  - 


0.0000 


0_q_c><^q^q^q_^_q^q_^__q-  Q  Q  g  Q-g-  Q--Q-  Q  Q-  Q  Q  Q  Q  Q 


Homogenization,  p  =  5.0(10) 


- v - 


O .  Conjugate  gradient,  pd  =  5.0(1 0)-4 

-▼ - Homogenization,  pd  =  1.2(1 0)-4 

Conjugate  gradient,  p  =  1 .2(1 0)-4 
Homogenization,  pd  =  5.5(1 0)'5 
Conjugate  gradient,  pd  =  5.5(1 0)'5 
Homogenization,  pd  =  3.1  (10)'5 
Conjugate  gradient,  p  =  3.1  (1 0)'5 
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Elastic  Stiffness 


Remarks 

-Minimum  energy  configurations 
comparable 

-Energy  compares  with  experimental 
result  of  3.6  eV  [Schultz,  1991] 

-Stiffness  affected  little  by  low  defect 
concentration 

-Strong  deformation-induced  anisotropy 
could  be  artifact  of  EAM  potential 


Zener  Anisotropy 


Screw  dislocation  configuration 


[1 111(11 0)  screw  dislocation  in  BCC  tungsten 
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unit  cell,  32256  atoms 


-screw  dislocation  core 
-relative  (111)  displacements  to 
initial  lattice  with  defect 
-1 1 1 -compression  (movie) 


-screw  dislocation  core 

-relative  (111 ^displacements  to  perfect 

BCC  lattice 

-null  applied  displacements 


-Three-fold  symmetry  of  W  dislocation  core  agrees  with  lattice  statics  solution  of  Vitek  [1976, 
Proc.  R.  Soc.  Lond.  A] 

-Tension-compression  asymmetry  of  stress  &  stiffness  agrees  with  experiments 

-Results  of  homogenization  method  cross-verified  using  sequential  lattice  statics  (Paradyn  code, 
S.  Plimpton) 
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Undeformed  mechanical  properties: 
W  unit  cells  with  periodic  dislocations 


m 
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Isolated  screw 
(lattice  curvature) 


Screw  dipole 

(counteracting  Burgers  vectors) 
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Elastic  constants:  dislocation  dipole 


0.00  0.05  0.10  0.15  0.20 


Pd  [1/nm2] 


Pd  [1/nm2] 


Screw  dislocation  with  elastic  shear 


28800  atom  simulation 
BCC  Tungsten 
[111]  Screw  dislocation 
1%  elastic  shear  strain  imposed 


Scaled  asymptotic  correction  v<k> 


g*  [eV/atom] 
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Screw  dislocation  -  [111]  tension 


Homogenization,  //  =  2.6(1 0)'2/nm2 
Conjugate  gradient,  p  =  2.6(1 0)'2/nm2 
Homogenization,  p  =  1.1(10)'2/nm2 
Conjugate  gradient,  p  =  1.1(1 0)'2/nm* 
Homogenization,  p  =  6.4(1 0)'3/nm2 
Conjugate  gradient,  =  6.4(1 0)'3/nm2 


0.005 
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Defect  Energy 


0.010  0.015 

/VI 


0.020  0.025 


-Comparable  minimum  energy  configurations 
attained  with  different  methods 
-Positive  coupling  of  defect  energy 
and  stretch 

-Stiffness  and  anisotropy  affected 
by  dislocations  ,  agrees  with  experimental 
trends  for  Cu  [Smith,  1953,  Phil.  Mag.  A] 
-Large  defect  densities:  ~10  nm  spacing 


Zener  Anisotropy 


[eV/ atom] 
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0.2-radian  twist  boundary:  [111]-tension 


Homogenization,/?  =  1.1(10)  /nm 

■O .  Conjugate  gradient, //  =  1.1(1 0)'2/nm2 

- ▼ - Homogenization,  p  =  8.2(1 0)'3/nm2 

- v -  Conjugate  gradient,  p  =  8.2(1 0)'3/nm2 


-►X 

(112) 
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-Comparable  minimum  energy  configurations 
attained  with  homogenization  method  and 
with  conjugate  gradient  technique 

-Negative  coupling  of  defect  energy 
and  applied  stretch 

-Stiffness  and  anisotropy  affected  significantly 
by  twist  boundary 
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Fixed  screw  dislocation  (atoms) 
superimposed  [110]-glide  (continuum) 


'Vv 


Defect  Configuration  Elastic-plastic  kinematics,  pure  shear: 
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Remarks: 

-Peierls-like  oscillations  depend  on  size  aA  of  unit  cell 
(controlled  by  distance  traveled  by  dislocation) 

-Elastic  stiffness  and  energy  changes  from  atomistics 
reflected  in  plasticity  simulation 


Plasticity  kinetics: 
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Resolved  Shear  Stress 


Elastic  +  Defect  Energy 


Elastic  Stiffness 


Summary 


•  Key  aspects  of  multiscale  theory 

-  Finite  deformation  with  evolving  intermediate  configuration 

-  Born  correction  (atoms)  linked  to  fine  scale  asymptotic  solution  (continuum) 

-  Separate  but  co-linear  coordinate  systems  at  each  scale 

-  No  need  to  match  boundary  interfaces  as  in  some  other  multiscale  methods 

-  Presently  restricted  to  periodic  defect  distributions  at  each  macro-element 

•  Implementation  and  examples 

-  Accurate  alternative  to  conjugate  gradient  minimization 

-  Convergence  and  performance  properties  not  fully  explored 

-  Demonstrated  effects  of  large  defect  concentrations  on  moduli  &  energy  of  W 

-  Static  formulation,  unable  to  predict  plasticity  kinetics  from  atomic  scale,  but 
should  be  possible  in  future  with  dynamic  formulation 

-  Incorporates  plasticity  kinematics  and  continuum-type  slip 

-  Mixed  domains  like  quasi-continuum  [Tadmor,  Ortiz,  Shenoy  et  al.]  possible 


